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Wingbeat  Shape  Modulation  for  Flapping-Wing 
Micro- Air- Vehicle  Control  During  Hover 

David  B.  Doman,*  Michael  W.  Oppcnhcimer,’  and  David  O.  SigthorssorC 
U.S.  Air  Force  Research  Laboratory,  Wright-Patterson  Air  Force  Base,  Ohio  45433-7531 

DOI:  10.2514/1.47146 

A  new  method  of  controlling  a  flapping- wing  micro  air  vehicle  by  varying  the  velocity  profiles  of  the  wing  strokes  is 
presented  in  this  manuscript.  An  exhaustive  theoretical  analysis  along  with  simulation  results  show  that  this  new 
method,  called  split-cycle  constant-period  frequency  modulation,  is  capable  of  providing  independent  control  over 
vertical  and  horizontal  body  forces  as  well  as  rolling  and  yawing  moments  using  only  two  physical  actuators,  whose 
oscillatory  motion  is  defined  by  four  parameters.  An  actuated  bob-weight  is  introduced  to  enable  independent 
control  of  pitching  moment.  A  general  method  for  deriving  sensitivities  of  cycle-averaged  forces  and  moments  to 
changes  in  wingbeat  kinematic  parameters  is  provided,  followed  by  an  analytical  treatment  for  a  case  where  the  angle 
of  attack  of  each  wing  is  passively  regulated  and  the  motion  of  the  wing  spar  in  the  stroke  plane  is  driven  by  a  split- 
cycle  waveform.  These  sensitivities  are  used  in  the  formulation  of  a  cycle-averaged  control  law  that  successfully 
stabilizes  and  controls  two  different  simulation  models  of  the  aircraft.  One  simulation  model  is  driven  by 
instantaneous  aerodynamic  forces  derived  from  blade-element  theory,  while  the  other  is  driven  by  an  empirical 
representation  of  an  unsteady  aerodynamic  model  that  was  derived  from  experiments. 


I.  Introduction 

EIGHT  and  volume  constraints  on  flapping-wing  micro  air 
vehicles  (MAV s)  drive  the  desire  to  use  as  few  actuators  as 
possible  while  achieving  a  high  level  of  controlled  maneuverability 
using  the  simplest  possible  controller.  The  best  way  of  achieving  this 
objective  is  an  open  question  and  represents  a  multidisciplinary 
design  problem  wherein  vehicle  controllability  must  be  considered 
from  the  outset.  Additionally,  the  issue  of  determining  an  adequate 
level  of  fidelity  for  model-based  control  design  for  flapping-wing 
aircraft,  which  are  influenced  by  unsteady  aerodynamic  phenomena, 
has  not  been  settled.  This  exposition  addresses  each  of  these  open 
problems. 

To  address  the  problem  of  reducing  the  number  of  required 
actuators  for  highly  manueverable  flight,  a  new  control  strategy, 
called  split-cycle  constant-period  frequency  modulation,  is  pre¬ 
sented,  which  shapes  each  wingbeat  velocity  profile  and  enables 
independent  control  of  vertical  and  horizontal  body  forces  and 
rolling  and  yawing  moments  using  only  two  physical  actuators. 
Thus,  four  degrees  of  freedom  can  be  controlled  using  two  physical 
actuators  whose  oscillatory  behavior  is  parameterized  by  four  inputs. 
The  overall  objective  of  achieving  a  high  level  of  maneuverability 
using  a  reduced  set  of  actuators  is  thereby  achieved,  which  represents 
an  advance  over  the  state  of  the  art  where  the  number  of  physical 
actuators  was  required  to  be  greater  than  or  equal  to  the  number  of 
independent  degrees  of  freedom  to  be  controlled  [1-7]. 

An  aircraft  concept  is  presented  that  is  equipped  with  three 
bimorph  piezoelectric  actuators.  Two  of  these  actuators  manipulate 
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the  wings  to  follow  split-cycle  wing  position  profiles,  while  the  third 
manipulates  a  bob-weight  that  allows  the  vehicle  center  of  mass  to 
shift  for  the  purpose  of  enabling  pitch  control.  The  result  is  the  ability 
to  independently  control  five  degrees  of  freedom  using  three  physical 
actuators.  Unlike  some  flapping-wing  MAV  designs  [5,6],  this 
aircraft  is  designed  to  be  capable  of  hovering  and  maneuvering  using 
only  two  wings  and  a  bob-weight,  without  the  need  for  conventional 
empennage  control  effectors. 

A  general  framework  is  presented  for  computing  a  control- 
oriented  model  for  flapping-wing  MAVs.  This  framework  is  used  to 
address  the  issue  of  determining  an  adequate  level  of  modeling 
fidelity  for  a  model-based  control  law  that  must  stabilize  and  control 
the  proposed  vehicle.  The  method  begins  by  constructing  a  cycle- 
averaged  model  of  the  aerodynamic  forces  and  moments  that  is 
similar  to  those  found  in  the  literature  [2,3, 8, 9],  However,  in  the 
present  method,  the  sensitivities  of  these  cycle-averaged  forces  and 
moments  with  respect  to  the  wingbeat  kinematic  parameters  that  can 
be  manipulated  are  computed  in  order  to  form  a  matrix  of  cycle- 
averaged  control  derivatives  that  can  be  used  for  control  synthesis. 
The  method  is  general  in  the  sense  that  the  aerodynamic  control 
derivatives  can  be  extracted  from  high-fidelity  numerical  aero¬ 
dynamic  models,  experimental  measurements,  or  simple  blade- 
element  models  depending  upon  the  desired  level  of  controller 
performance.  In  the  interest  of  determining  the  level  of  adequate 
fidelity  for  a  model-based  controller,  an  analytical  blade-element 
aerodynamic  model  is  presented  for  the  proposed  aircraft.  The 
control  derivatives  are  then  analytically  derived  and  used  to  design  a 
controller.  These  sensitivities  are  used  in  the  formulation  of  a  cycle- 
averaged  control  law  that  successfully  stabilizes  and  controls  two 
different  simulation  models  of  the  aircraft.  One  simulation  model  is 
driven  by  instantaneous  aerodynamic  forces  derived  from  blade- 
element  theory,  while  the  other  is  driven  by  an  empirical  repre¬ 
sentation  of  an  unsteady  aerodynamic  model  that  was  derived  from 
experiments.  Simulation  results  indicate  that  a  cycle-averaged  blade- 
element-model-based  controller  is  robust  to  the  aerodynamic  mod¬ 
eling  error  introduced  by  the  control-oriented  modeling  assumptions. 
Additionally,  it  is  noted  that  stability  and  tracking  performance  for 
flapping- wing  aircraft  must  be  interpreted  in  terms  of  how  the  aircraft 
is  targeted  to  equilibrium  points  that  it  will  ultimately  orbit.  Limit 
cycles  about  target  set  points  are  unavoidable  [10-12]  because 
flapping  wings,  by  their  very  nature,  produce  periodic  forces  and 
moments. 

Finally,  by  computing  an  analytical  control-oriented  model  from 
blade-element  theory,  a  direct  link  between  the  aircraft  design  and  the 
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vehicle’s  controllability  can  be  established.  This  feature  allows 
designers  to  consider  controllability  early  in  the  design  process  and 
enables  insights  into  control  strategies  that  might  be  obfuscated  by 
purely  numerical  representations  of  the  aerodynamics.  Insights  into 
vehicle  physics  and  design  are  highlighted  throughout  the  analytical 
development. 

II.  Vehicle  Concept 

The  first  takeoff  of  an  insect-scale  flapping-wing  MAV  was 
achieved  by  an  aircraft  called  RoboFly  that  was  developed  at  Harvard 
University  by  Wood  [13].  A  key  feature  that  led  to  its  successful  first 
flight  is  that  the  vehicle  was  minimally  actuated.  RoboFly  uses  a 
single  bimorph  piezoelectric  actuator  to  impart  symmetric  motion  to 
two  wings.  The  wing  angle  of  inclination  with  respect  to  the  stroke 
plane  is  regulated  by  passive  wing  rotation  joints  resulting  in 
planform  motion  that  is  similar  to  that  of  a  dipterian  insect.  This 
passive  rotation  of  the  planform  about  the  spar  eliminates  the  need  for 
added  weight  and  complexity  associated  with  actuators  that  actively 
rotate  the  planform  about  the  spar.  The  first  flight  at  Harvard  resulted 
in  unregulated  flight  up  a  pair  of  wires  that  constrained  the  vehicle 
motion  to  vertical  translation.  In  the  present  paper,  a  vehicle  concept 
similar  to  the  Harvard  Robofly  is  proposed,  along  with  a  control 
strategy,  that  would  enable  controlled  six-degree-of-freedom 
(DOF)  flight  of  the  fuselage  without  the  need  for  artificial  motion 
constraints. 

A.  Design  Features 

A  diagram  of  the  proposed  vehicle  is  shown  in  Fig.  1 .  The  main 
differences  between  the  proposed  vehicle  and  the  RoboFly  are 
that  the  new  control-configured  aircraft  would  be  equipped  with 
independently  actuated  wings  and  an  actuated  bob-weight  that 
allows  the  vehicle  center  of  gravity  to  be  manipulated  for  control 
purposes.  The  wings  rotate  about  hinges  at  the  wing  root.  A  linkage- 
based  transmission  translates  the  tangential  motion  of  the  tip  of  a 
bimorph  piezoelectric  actuator  into  rotational  motion  of  the  wings 
in  the  stroke  plane.  The  linkage  elements  are  designed  to  achieve 
impedance  matching  between  the  wing  and  actuator  forces  as  well  as 
amplify  the  relatively  small  motion  of  the  tip  of  the  bimorph  strip  into 
large  angular  displacements  of  the  wing  root.  The  planforms  are 
connected  to  the  movable  wing  roots  by  a  limited  hinge  or  flexure 
joint  that  provides  for  passive  rotation  of  the  wing.  This  hinge  allows 


the  wing  to  passively  flip  over  as  the  wing  reverses  direction  at  the 
end  of  each  stroke  and  enables  the  chord  to  rotate  about  the  axis  of  the 
spar  in  a  manner  that  approximates  the  wing-twisting  motion  that  has 
been  observed  in  dipterian  insects  [14],  As  the  wing  rotates  through 
the  stroke  plane,  dynamic  pressure  acting  on  the  wing  tends  to  cause 
it  to  feather  into  the  wind;  however,  as  shown  in  Fig.  2,  a  flexure  hinge 
joint  is  designed  such  that  the  spar  and  root  structure  interfere  at  an 
angle  set  by  the  designer  that  prevents  the  wing  from  overrotating. 
This  interference  causes  the  wing  to  hold  a  constant  angle  of  attack 
relative  to  the  stroke  plane  once  a  critical  dynamic  pressure  is 
reached.  The  bimorph  piezoelectic  actuators  and  the  carbon  fiber 
substrate  to  which  they  are  mounted  are  cantilevered  to  the  fuselage. 
In  the  case  of  the  wings,  the  actuator  assembly  and  linkage,  together 
with  the  wing,  form  a  spring-mass-damper  system  that  has  a  known 
resonant  frequency.  The  actuator  assembly  for  the  bob-weight  is 
simply  used  to  control  the  mean  pitch  attitude  of  the  fuselage.  In  the 
single  actuator  Harvard  experiments,  the  dynamic  system  was  driven 
at  resonance  [15]  for  maximum  energy  efficiency  to  achieve  flight.  In 
the  present  problem,  each  wing-linkage-actuator  system  is  nominally 
driven  at  a  hover  frequency  that  is  defined  as  the  frequency  at  which 
the  cycle-averaged  lift  is  equal  to  the  aircraft  weight.  The  term  hover 
is  to  be  inteipreted  in  a  time-averaged  sense  because,  as  we  will  show, 
the  aircraft  is  in  constant  motion  due  to  the  periodic  nature  of  the 
forces  and  moments  produced  by  the  flapping  wings.  The  closest 
approximation  of  hover  that  can  be  achieved  is  that  of  a  high- 
frequency,  low-amplitude  limit  cycle  about  a  mean  position  [12], 


Downstroke  End  of  Stroke  Upstroke 


Fig.  1  General  assembly  of  a  control-configured  MAV  equipped  with  independently  actuated  wings  and  an  actuated  bob-weight. 
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B.  Control  Strategy 

The  coordinate  frames  and  wing  angles  used  in  the  development  of 
the  control  strategy  are  illustrated  in  Fig.  3.  The  vehicle  is  designed  to 
mimic  dipterian  insect  flight;  thus,  the  angular  velocity  of  each  wing 
in  its  respective  stroke  plane  <f>  must  be  a  time- varying  function  that  is 
equal  to  zero  at  the  extreme  limits  of  wing  position.  It  is  assumed 
that  the  wing  position,  in  the  stroke  plane,  can  be  controlled  directly 
via  a  high  bandwidth  piezoelectric  bimorph  actuator.  To  begin  the 
development,  a  simple  cosine  function  is  selected  to  drive  the  stroke 
plane  rotation  angle  of  each  wing: 

< p  =  cos  cot  (1) 

The  maximum  value  of  this  angle  occurs  when  either  wing  reaches  its 
maximum  displacement  towards  the  ventral  side  of  the  vehicle. 
Assuming  that  the  frequency  of  the  oscillating  wing  is  held  constant 
over  the  segment  of  interest  for  each  wingbeat  cycle,  the  angular 
velocity  of  the  wing  is  given  by 

cf>  =  —co  sin  cot  (2) 

Note  that  the  units  of  <f>  are  in  radians  and  that  the  amplitude  of  the 
wing  rotation  in  the  stroke  plane  is  taken  to  be  ±  1  rad,  which  closely 
approximates  the  rotation  limits  used  in  the  Harvard  RoboFly  and 
those  observed  in  some  flying  insects  such  as  the  hawkmoth 
Manduca  sexta  [16].  The  frequency  of  the  oscillator  that  drives  the 
actuator  is  selected  as  one  of  the  control  input  variables  that  is  used 
to  enable  multidegree-of-freedom  flight;  however,  this  frequency  is 
allowed  to  change  at  key  points  in  the  wingbeat  cycle,  namely  at 
the  ends  of  the  stroke  when  <f>  =  1,0  =  0  and  0  =  — 1,0  =  0.  If  the 
wingbeat  frequency  is  held  constant  over  each  cycle  and  the  angles  of 
attack  for  each  stroke  are  equivalent,  then  the  cycle-averaged 
longitudinal  force  along  the  x  body  axis  is  always  finite  and  positive, 
while  the  cycle-averaged  z  body  force  is  always  zero.  Such  con¬ 
ditions  make  it  possible  to  achieve  a  practical  condition  that  approxi¬ 
mates  hover  and  allows  altitude  tracking  [12];  however,  to  enable  a 


vehicle  to  achieve  multidegree-of-freedom  flight,  an  additional 
control  input  variable  must  be  introduced. 

To  achieve  a  nonzero  cycle-averaged  body  force  in  the  z  direction, 
using  wings  equipped  with  passive  wing  rotation  joints  as  described 
above,  the  stroke  plane  angular  velocity  of  each  wing  must  be 
temporally  asymmetric  over  the  wingbeat  cycle;  i.e.,  the  wing 
velocity  profiles  of  the  upstroke  and  downstrokes  must  be  different. 
Ennos  [17]  presents  biological  evidence  for  differing  upstroke  and 
downstroke  velocities.  One  way  to  induce  temporal  asymmetry 
across  each  wingbeat  is  to  parameterize  the  waveforms  that  define  the 
motion  of  each  wing  in  the  stroke  plane  in  tenns  of  temporally 
symmetric  and  asymmetric  components.  The  new  technique  that  is 
proposed  to  produce  this  desired  behavior  is  referred  to  as  split-cycle 
constant-period  frequency  modulation.  The  fundamental  idea  is  to 
piece  together  two  cosine  waves  of  different  frequencies  over  one 
complete  wingbeat  cycle,  i.e.,  an  upstroke  and  a  downstroke.  The 
period  of  the  cycle  is  fixed  by  the  symmetric  or  fundamental 
frequency,  while  the  speed  of  each  stroke  within  the  cycle  is  governed 
by  an  asymmetric  frequency  or  split-cycle  parameter. 

A  detailed  discussion  of  the  mathematical  construction  of  these 
split-cycle  waveforms  is  provided  by  Doman  et  al.  [18]  and 
Oppenheimer  et  al.  [19].  In  short,  the  split-cycle  control  strategy 
generates  time-varying  upstroke  and  downstroke  wing  position 
commands  that  are  defined  by 

0t/(O  =  cos[(oj  —  5)r],  0  <t<- — — —  (3) 


it  2t r 

0d(O  =  cos[(w  +  a)t  +  £],  - - K-f< —  (4) 

(co  —  o)  co 

where 


and 


Fig.  3  Relationship  between  body,  root,  spar,  left-,  and  right-wing-planform  axis  systems  on  upstroke. 
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S>  0 


Normalized  Time 


Fig.  4  Split-cycle  waveforms. 


htenrelizecJTirre 


6<0 


—2jr8 

co-28 


(6) 


GW))  & 


(9) 


Thus,  the  waveform  that  drives  the  angular  position  of  each  wing  in 
the  stroke  plane  is  completely  defined  by  the  symmetric  frequency  co 
and  the  asymetric  frequency  or  split-cycle  parameter  8,  where  <5  <  |. 
The  stroke  plane  angular  velocity  of  each  wing  over  each  split-cycle 
interval  is  given  by 

4>  u(t)  =  -(co  -  8)  sin[(w  -  3)/],  0<f<  n  (7) 

((o-o) 

<Pd({)  =  -(w  +  ff)sin[(w  +  ff)r+f].  n  <  t  <  —  (8) 

(co  —  8)  co 


Since  the  introduction  of  the  split-cycle  parameter  changes  the 
frequency  of  the  cosine  wave  when  (co  —  8)t  =  n  in  each  cycle,  it  is 
convenient  to  split  the  integral  as  follows: 


r-M 
G  ~  lit 


fmG(cf>u(t))dt+  f“  G(Mt))to 


(10) 


The  control  derivatives  associated  with  a  generalized  cycle- 
averaged  force  are  calculated  by  computing  the  partial  derivative 
with  respect  to  each  control  input  variable  about  the  hover  frequency 
coa.  The  control  derivatives  at  hover  are 


Figure  4  shows  the  wingbeat  position  profile  for  8  >  0  and  8  <  0 
along  with  a  symmetric  wingbeat.  The  split-cycle  parameter  8  for 
each  wing,  along  with  co  for  each  wing,  become  control  parameters. 
Hence,  the  waveforms  that  drive  the  left  and  right  wing  positions  are 
defined  by  the  following  control  parameters:  coLW,  8LW,  coRW,  and 
8rw,  respectively. 

The  control  strategies  proposed  in  this  manuscript  are  based  on  the 
assumption  that  the  bandwidth  of  the  closed-loop  fuselage  control 
system  is  much  less  than  the  flapping  frequency  required  to  achieve  a 
quasi-hover  condition.  If  nonoscillatory  control  forces  and  inertial 
measurements  were  available  for  feedback,  one  could  use  con¬ 
ventional  control  design  techniques  to  synthesize  feedback  control 
laws  that  produce  favorable  closed-loop  responses.  However, 
because  of  the  nature  of  the  aerodynamic  forces  produced  by  flapping 
wings,  time- varying  high-frequency  oscillatory  forces  and  moments 
must  be  manipulated  to  impart  desired  motion  to  the  fuselage.  Under 
this  constraint,  the  tracking  behavior  of  the  fuselage  will,  under  the 
best  of  conditions,  be  characterized  by  low-amplitude  and  high- 
frequency  limit  cycles  about  desired  set  points.  Towards  this  end,  the 
relationships  between  the  cycle-averaged  forces  and  moments  and 
the  control  input  parameters  that  define  the  characteristics  of  the 
oscillators  that  drive  the  wing  positions  are  computed.  Feedback 
control  laws  based  on  cycle-averaged  forces  and  moments  will  be 
designed  to  allow  a  vehicle  to  track  desired  angular  and  spatial 
positions  in  a  mean  sense. 

The  general  method  used  to  compute  the  cycle-averaged 
aerodynamic  forces  and  their  sensitivities  to  variations  in  the  param¬ 
eters  that  define  the  wingbeat  kinematics  will  now  be  presented. 
The  method  is  general  in  the  sense  that  the  aerodynamic  control 
derivatives  can  be  extracted  from  high-fidelity  numerical  aero¬ 
dynamic  models,  experimental  measurements,  or  simple  blade- 
element  models  depending  upon  the  desired  level  of  controller 
performance.  G(t)  is  defined  to  be  a  generalized  force  that  is  aligned 
with  a  principal  body-axis  direction.  Thus,  G(t)  may  represent 
either  a  force  or  a  moment.  To  compute  the  cycle-averaged  gener¬ 
alized  force  associated  with  each  wing,  an  integral  of  the  following 
form  is  evaluated: 


JG_ 

3  «i.w 


(ID 


wLW=wo.^LW=0 


Gs,„,  = 


a  g 
38,  w 


(12) 


^LW=(t)o^LW=0 


3  G 


(13) 


(Orw=<Do,<5rw=0 


G 


<5rw 


dG 

3SRW 


coRW=(oo,SRW=0 


where  the  hover  frequency  [12]  is  given  by 


2  mg 
pIACL(a) 


(14) 


(15) 


These  control  derivatives  can  be  used  to  design  the  control 
allocation  portion  of  a  six-DOF  control  law  for  the  proposed 
flapping-wing  MAV.  The  above  calculations  could  be  carried  out 
numerically  if  high-fidelity  unsteady  computational  fluid  dynamic 
simulations  or  experimental  measurements  were  available  and  it 
was  determined  that  lower-fidelity  methods  produced  model-based 
controllers  that  were  incapable  of  providing  adequate  performance. 
As  will  be  shown  throughout  the  remainder  of  the  exposition,  a 
blade-element-based  aerodynamic  model  is  sufficiently  accurate  to 
produce  controllers  that  provide  adequate  performance  when 
applied  to  a  simulation  model  that  is  driven  by  unsteady  aero¬ 
dynamics.  The  blade-element  model  has  the  added  advantage  that 
cycle-averaging  calculations  can  be  carried  out  analytically  and  the 
resulting  model  shows  a  clear  and  direct  tie  between  key  vehicle 
design  parameters  and  its  controllability. 
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III.  Instantaneous  Blade-Element  Model 

Strictly  speaking,  all  mathematical  models  of  physical  systems  are 
idealized  representations  of  systems  that  cannot  physically  exist.  The 
critical  issue  for  the  control  engineer  is  to  formulate  a  model  of  a 
physical  system  of  interest  that  is  of  suitable  fidelity,  which  enables  a 
control  system  to  be  designed  that  is  robust  to  errors  resulting  from 
the  modeling  process.  A  control-oriented  model  should  be  compu¬ 
tationally  efficient  so  that  dynamic  simulations  of  the  closed-loop 
physical  system  run  in  a  reasonable  period  of  time;  furthermore,  it  is 
desirable  for  such  models  to  provide  a  direct  link  between  vehicle 
controllability  and  the  aircraft  design  parameters. 

The  purpose  of  this  development  is  to  derive  an  idealized 
analytical  model  that  will  allow  one  to  study  the  dynamic  behavior  of 
the  aircraft  under  consideration  and  to  formulate  and  ultimately  test 
the  proposed  control  strategies  using  nonreal-time  simulations  that 
relax  many  of  the  key  assumptions  used  to  derive  the  analytical 
model.  The  following  assumptions  are  made  in  the  formulation  of 
this  analytical  control-oriented  model:  no  aerodynamic  interactions 
exist  between  the  left  and  right  wings,  nor  between  the  wings  and  the 
fuselage;  two-dimensional  sectional  aerodynamic  coefficients  are 
known  and  constant  throughout  each  stroke;  passive  wing  rotation 
joints  are  on  a  limit  when  the  corresponding  wing’s  stroke  plane 
angular  velocity  is  nonzero;  aerodynamic  forces  and  moments  are  the 
sole  result  of  wing  motion,  i.e.,  only  flight  conditions  in  the 
neighborhood  of  hover  are  considered;  and  piezoelectric  actuator 
bandwidth  is  much  greater  than  the  wingbeat  frequency  required  to 
hover. 

The  dynamic  analysis  of  the  vehicle  requires  the  use  of  several 
coordinate  frames.  This  is  because  the  signs  of  aerodynamic  forces 
and  moments  differ  from  the  left  wing  to  the  right  wing.  Furthermore, 
the  signs  of  the  forces  and  moments  are  dependent  upon  the  direction 
of  the  wing  stroke,  i.e.,  upstroke  or  downstroke.  The  aerodynamic 
forces  and  moments  can  be  conveniently  written  in  certain  inter¬ 
mediate  frames;  however,  it  will  be  necessary  to  transform  them  into 
a  body-fixed  and  ultimately  an  inertial  coordinate  frame  in  order  to 
write  the  equations  of  motion  of  the  fuselage.  The  coordinate  frame 
and  angle  definitions  are  illustrated  in  Fig.  3.  A  comprehensive  table 
containing  the  definitions  of  the  numerous  coordinate  frames  and 
rotation  matrices  is  provided  by  Doman  et  al.  [18]. 


A.  Instantaneous  Aerodynamic  Forces  and  Centers  of  Pressure  in 
Wing  Planform  Frames 

The  aerodynamic  forces  are  derived,  using  blade-element  theory, 
for  wings  that  have  two  degrees  of  freedom,  namely,  angular 
displacement  0(  t)  about  the  wing  root  in  the  stroke  plane  and  angular 
displacement  of  the  planform  about  a  passive  rotation  hinge  joint, 
which  is  equivalent  to  wing  angle  of  attack  a ,  in  still  air.  The  planform 
is  taken  to  be  a  rigid  flat  plate  and  the  lift  and  drag,  produced  by  each 
wing  [18],  can  be  expressed  as 

L  =  kL(p(t)2  D  =  kD<p(t)2  (16) 

where 

kL=^CL(a)lA  kD  =  ^CD(a)IA  (17) 

and  IA  is  the  area  moment  of  inertia  of  the  planform  about  the  axis  of 
the  root-hinge,  p  is  the  atmospheric  density,  and  CD(a)  and  CL(a) 
are  the  lift  and  drag  coefficients.  Empirical  expressions  for  lift  and 
drag  coefficients  [20],  which  were  obtained  from  low-Reynolds- 
number  experiments,  are  used  in  this  work: 

CD(a)=  1.92  —  1.55  cos(2. 04a  —  9.82) 

CL(a)  =  0.225  +  1.58  sin(2.13o!-  7.2)  (18) 


B.  Expression  of  Aerodynamic  Forces  and  Centers  of  Pressure 
in  the  Body  Frame 

Using  the  relationships  between  the  body,  roots,  spars,  upstroke- 
planform,  and  downstroke-planform  axis  systems  defined  in  Doman 
et  al.  [18],  the  instantaneous  values  of  lift  and  drag  on  each  wing  can 
be  transformed  into  the  body-axis  coordinate  frame.  If  the  air  mass  is 
quiescent,  then  the  relative  wind  is  parallel  to  the  stroke  plane,  which 
is  coincident  with  the  x-y  planes  of  the  local  spar  frames.  Therefore, 
the  lift  and  drag  forces  are  conveniently  expressed  in  the  spar-carried 
coordinate  frames,  which  can  be  transformed  to  the  body  frame.  The 
aerodynamic  forces  associated  with  each  wing  and  stroke  expressed 
in  both  the  spar  and  body  frames  are  summarized  in  Table  1 . 

The  center  of  pressure  of  a  wing  rotating  about  a  fixed  hinge  line 
can  be  computed  using  blade-element  theory  [18].  The  center  of 
pressure  of  each  wing  can  be  conveniently  expressed  in  the  local 
wing  planform  frame  associated  with  each  stroke.  To  calculate  the 
instantaneous  center  of  pressure  in  the  body  frame,  three  coordinate 
frame  rotations  and  one  translation  must  be  performed.  The  resulting 
expressions  for  the  centers  of  pressure  associated  with  each  wing  and 
stroke,  expressed  in  the  body  frame,  are  summarized  in  Table  2.  In 
these  expressions,  and  y^p  are  the  x  and  v  locations  of  the  center 
of  pressure  in  the  local  wing  planform  frame,  Arf  is  defined  as  the 
position  vector  from  the  origin  of  the  body  frame  to  the  origin  of  the 
right  wing  root  coordinate  system,  i.e.. 


\rB  —  rB 

~  roRWR 

Ar*  A  rB 

ArL  —  r°LWR 


a  4  f  AzJ 

AxB  -f  Azf 


(19) 


and  the  width  of  the  vehicle  is  defined  as  w. 


Table  1  Aerodynamic  forces  expressed  in  local  spar  and  body  frames 


Force 

Local  spar  frame 

Body  frame 

RW  upstroke 

|7  RWS  _ 

1  RWU  ~~ 

£*RWU 

0 

_  — ^RWU  _ 

F  B  — 

r  RWU  — 

7-rwu 

— Drwu  sin  0RWU  ( t ) 

.  DRwucos<()RWU(f)  _ 

RW 

downstroke 

f?  RWS  _ 

1  RWD  — 

— ^RWD 

0 

_  ^RWD _ 

g» 

II 

7-rwd 

Dr wd  sin0RWD(f) 

— Drwd  cos  0RWD  ( t) . 

LW  upstroke 

*1 

rr 

5=5- 

CjC/5 

II 

— ^LWU 

0 

_  — ^LWU  _ 

T 

a 

II 

7-lwu 

D  lwu  sin  0lWU(7  ) 

_  Dlwu  cos  <A.wu(f)  _ 

LW 

downstroke 

<< 

O'* 

II 

f\wD 

0 

 ^LWD  

o 

II 

^LWD 

f^LWD  sin  0LWD  W 
^LWD  cos  0LWd(O  _ 

Table  2  Centers  of  pressure  expressed  in  body  frame 


CP  location 

Body  frame  expression 

RW  upstroke 

r B  — 

1  cpRWU  — 

x)Jp  sin  a  +  Axf 
sin  <Pr w  cos  “  +  cos  0RW  +  f 

-a*p  cos  < pRW  cos  a  +  y™p  sin  0RW  +  Azf  _ 

RW 

downstroke 

r  B  — 

r  cpRWD  — 

a™  sin  a  +  Axf 

-.v*p  sin  <pRW  cos  a  +  >AP  cos  0RW  +  f 
_  cos  (/> RW  cos  a  +  y™  sin  0RW  +  Azf  _ 

LW  upstroke 

rB  _ 

1  cpLWU  — 

sin  a  +  A.rf 

-x*p  sin  tpLW  cos  a  -  y*p  cos  </>LW  -  f 
_  cos  0LW  cos  a  +  sin  0LW  +  Azf  _ 

LW 

downstroke 

rB  _ 

1  cpLWD  “ 

x™  sin  a  +  A.rf 
sin  <pLW  cos  a  -  y™  cos  0LW  -  f 
 xWp  cos  0LW  cos  a  +  y™  sin  0LW  +  Azf  

where  a  is  in  degrees. 
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C.  Aerodynamic  Moments  in  Body  Frame 

The  expressions  for  the  aerodynamic  moments  associated  with 
each  wing  and  stroke  are  given  by 


pin 

«Si 


sin2[(&>  +  a)t  +  £]  dr  = 


2(&>  +  o') 


(28) 


M B  —  rB  x  F5 
RWU  —  1  cpRWU  x  r  RWU 

Mlwu  =  rfpLwu  x  Flwu 


Mfi  —  rB  x  Ffl 
ivlRWD  —  1  cpRWD  A  rRWD 

MfwD  =  rcpLWD  X  FlWD 


(20) 


Carrying  out  the  cross-product  operations  and  substituting  the  values 
from  Tables  1  and  2  into  the  expressions  in  Eq.  (20)  yields 


/In 

“  sin2[(&)  +  a)t  +  f]  sin(cos[(w  +  a)t  +  £])  dr  =  0  (29) 


/In 
> 


h=  I  "  sin2[(a)  +  a)f  +  ^]cos(cos[(n)  +  g)r  +  ^])dr=  ^  ^  (30) 

co  +  a 


M 


B  _ 
RWU  — 


^RWu[iVcpP  +  2  cos0RW  +  Aif  Sttt'/’Rw] 

{^rwu[v^P  sin0RW  +  Ati£]  —  ^RWuA-tf  cos0rw  ~  [^rwu  cos  a  +  DRWV  sina]x*p  cos  0RW} 
{— Atwu  A.vjj  sin  0rw  —  ^Rwulf  +  T^P  cos  ‘/'rw]  —  I^rwu  cos  a  +  ^rwu  sin  “]-v^P  sin  0rw} 


M 


B  _ 
RWD  — 


— ^rwd[v^P  +  f  cos^rw  +  A sin0RW] 

{^RWD[y^P  sin'/'Rw  +  Az£]  +  Z)rwdAv|cos0rw  +  [Z,RWDcosa  +  Drwd  sina]**15  cos  0RW} 
{^rwdAy|  sin(/)RW  —  TRWD[f  +  T^Pc°s0rw]  +  [Arwd  cos  a  +  Drwd  sina]v^p  sin(/>RW} 


(21) 


(22) 


M 


B  _ 
LWU  — 


Dlwu[  TcpP  —  f  cos0LW  —  Az£  sin0LW] 

{^lwu[.V^p  sin  0LW  +  Az£]  —  -Dlwu  Ax£  cos  0lw  —  [Llwu  cos  a  +  Z)LWU  sin  cos  0LW } 
{^lwuAy£  sin0LW  +  Li,wu[f  +  TcpP  cos  0Lw]  +  [^lwu  cos  a  +  Dlwu  sinaf]x^p  sin(/>LW} 


(23) 


M 


B  _ 
LWD  — 


^lwdL)'cpP  +  f  cos0LW  +  A Zi  sin0LW] 

{^lwdLvcpP  sin  0lw  +  Az£]  +  Dlwd  Ax£  cos  (f>LW  +  [Llwd  cos  a  +  Dlwd  sin  a]x^p  cos  0LW} 
{— ^lwdAy£  sin</>LW  +  f-LWDlf  +  TcpP  cos^>LW]  —  [Tlwd  cos  a  +  Dlwd  sin  a]x*p  sin0LW} 


(24) 


Equations  (21-24)  provide  the  expressions  for  the  instantaneous 
aerodynamic  moments  generated  by  each  wing  at  any  point  in  a 
wingbeat  cycle  in  the  body-axis  coordinate  frame. 

IY.  Cycle-Averaged  Blade-Element  Model 

To  calculate  the  cycle-averaged  forces  and  moments  produced  by 
wings  that  follow  a  split-cycle  position  prohle,  it  is  necessary  to 
evaluate  numerous  integrals.  Many  of  the  integrands  are  of  the  form 
cos(cos  cot )  or  sin(cos  cor)  and  have  no  indefinite  integral  solutions. 
These  terms  arise  as  a  result  of  coordinate  transformations  between 
local  and  body  coordinate  frames,  coupled  with  the  fact  that  the  time- 
varying  wing  position  is  defined  in  terms  of  cosine  wave  segments. 
Fortunately,  definite  integrals  involving  such  functions  exist  over  the 
intervals  of  interest  for  the  present  problem  and  can  be  derived  from 
results  presented  in  Gradshteyn  and  Ryzhik  [21].  The  solution  to 
many  of  these  definite  integrals  involves  a  Bessel  function  of  the  first 
kind,  ■/,(•),  and  a  Struve  function.  The  integrals  involving  Struve 
functions  evaluate  to  zero  over  the  limits  of  integration  that  arise  in 
the  formulation  of  this  particular  problem.  For  convenience,  the 
solutions  to  the  definite  integrals  that  are  required  to  compute  the 
cycle-averaged  forces  and  moments  are  provided  below: 

A  =  ‘  sin2[(«  -  5)/]  At  =  U  (25) 

Jo  2(a)  —  0) 

I2  =  f  sin2[(<w  —  5)r]  sin(cos[(o>  —  6)?])  dr  =  0  (26) 

Jo 


A.  Cycle-Averaged  Forces 

In  this  section,  expressions  for  cycle-averaged  aerodynamic  forces 
are  derived  and  expressed  in  terms  of  the  split-cycle  waveform 
control  input  parameters  that  govern  the  wing  kinematics  and  time 
invariant  constants.  The  instantaneous  forces  are  integrated  with 
respect  to  time  as  each  wing  operates  over  an  upstroke  and 
downstroke. 


1.  X  Force 

Substituting  the  expression  for  the  instantaneous  x  body  force 
from  Table  1  into  Eq.  (10)  produces 


y  B  _ 

^  RW  — 


frjRW 
2  n 


(31) 


Substituting  Eqs.  (7),  (8),  and  (16)  into  Eq.  (31)  and  noting  that  the 
time-varying  functions  under  the  integral  signs  are  of  the  form  /,  and 
/4  as  given  by  Eqs.  (25)  and  (28),  respectively,  Eq.  (3 1 )  can  be  written 
as 


A rw  —  ^Wrw  —  ^Rw)2A  +  (tt’RW  +  Orw)2A]  (32) 

or  simply 

^rw  =  L  ^RW  (2wrw  ~  <Vw  +  °rw)  (33) 

Following  a  similar  procedure  for  the  left  wing,  it  can  be  shown  that 


73  =  f 


sin2[(a>  —  5)/]  cos(cos[(&>  —  <5)r])  df  = 


nJ  i(l) 
(a>-8) 


(27) 


yS  _  kL«> LW  n 
A  LW  —  A  yZCOL 


■  <Vw  +  °Lw) 


(34) 
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Note  that  both  XRW  and  Xfw  are  positive  quantities.  At  a  hover 
condition,  where  the  x  body-axis  is  normal  to  the  surface  of  the  Earth, 
the  forces  produced  by  both  wings  act  to  counter  the  vehicle  weight. 


2.  Y  Force 

To  compute  the  cycle-averaged  y  body  force  for  each  wing,  the 
expression  for  the  instantaneous  y  body  aerodynamic  force  from 
Table  1  is  substituted  into  Eq.  ( 10)  to  obtain 


yfw  =  ^  [  -ZWO  sintfcwW]  dt 

/2ji 

DrwdO)  sin[0RW(r)]  d? 

“RW-^RW 


(35) 


Substituting  Eqs.  (3),  (4),  (7),  (8),  and  (16)  into  Eq.  (35)  and  noting 
that  the  time-varying  functions  under  the  integral  signs  are  of  the 
form  /2  and  I5  as  given  by  Eqs.  (26)  and  (29),  respectively,  Eq.  (35) 
becomes 


1.  Rolling  Moment 

The  cycle-averaged  x  body  moment  can  be  written  in  terms  of 
Eq.  (10)  as 


Mx  rw  : 


^RW 

2n 


I  "RW^RW 

.Jo 


AZ.vrwuW  dt  + 


/: 


Mr, 


Jt)dt 


"RW^RW 


(42) 


Substituting  Eqs.  (3),  (4),  (7),  (8),  (16),  (21),  and  (22)  into  Eq.  (42) 
and  noting  that  the  integrals  are  of  the  form  given  in  Eqs.  (25-30), 
Eq.  (42)  can  be  written  as 


3Zcrw  —  ■ 


2n 


-{ 


(“rw  —  ^rw)2  T^PZ|  +  Azf/2  +  —I3 


~  (<Wrw  +  Orw)‘ 


TcpP/4  +  Az|/5  +  ~grh 


(43) 


Substituting  the  results  for  the  definite  integrals  I{  —  /6  and 
simplifying  yields 


yfw  =  0  (36) 

Following  a  similar  procedure  for  the  left  wing,  it  can  be  shown  that 

y£w  =  0  (37) 

Physically,  this  result  arises  from  the  fact  that  over  a  single  stroke,  the 
component  of  drag,  pointing  in  the  y  body  direction,  points  in  the 
positive  direction  for  the  same  amount  of  time  as  it  does  in  the 
negative  direction  and  with  the  same  magnitude. 


KU  =  _kDu  rwOW  +  ^wp  +  w7i(1)]  (44) 

Following  a  similar  procedure  for  the  left  wing,  it  can  be  shown  that 
TTB  ^D^LwlAw  +  °Lw)  r  wp  ,  , 

Mx  LW  = - 4 - [Top  +  wJ  1  (1)]  (45) 

Note  that  without  split-cycle  frequency  modulation,  i.e., 
i5rw  =  crRW  =  £lw  =  aLW  =  0,  it  would  not  be  possible  to  generate 
nonzero  cycle-averaged  rolling  moments  on  this  vehicle. 


3.  Z  Force 

The  calculation  of  the  cycle-averaged  force  in  the  z  body  direction 
follows  a  similar  procedure.  Substituting  the  expression  for  the 
instantaneous  z  body  aerodynamic  force  from  Table  1  into  Eq.  (10) 
yields 


yB  _  ^Rw 

zrw  — 

2jr 


D 


RWU 


(f)cos[0RW(t)]d? 


/2  Tt 

~ Z)RWd(0  c°s[0Rw(t)]  dt 

'"RW^RW 


(38) 


Substituting  Eqs.  (3),  (4),  (7),  (8),  and  (16)  into  Eq.  (38)  and  noting 
that  the  time-varying  functions  under  the  integral  signs  are  of  the 
form  Z3  and  Z6  as  given  by  Eqs.  (27)  and  (30),  respectively,  Eq.  (38) 
can  be  written  as 


Zrw  —  °2^W  [(WRW  —  ^rw)“Z3  —  (wRW  +  ctrw)2/6]  (39) 


or  simply 


_  ~^D-Z[(l)ft)RW 


(Vw  +  °rw)  (40) 


Following  a  similar  procedure  for  the  left  wing,  it  can  be  shown  that 


Z 


B  _ 
LW  — 


—kpJ  |(l)ft>LW 

2 


OW  +  °Lw) 


(41) 


2.  Pitching  Moment 

The  instantaneous  y  body  moment  from  Eq.  (10)  is 


rw  — 


^RW 

2n 


f 

Jo 


MyR\fnj(t)  dt  + 


L 


M, 


yRWD 


(t)dt 


wRW-aRW 


(46) 


Substituting  Eqs.  (3),  (4),  (7),  (8),  (16),  (21),  and  (22)  into  Eq.  (46) 
and  noting  that  the  integrals  are  of  the  form  given  in  Eqs.  (25-30), 
Eq.  (46)  can  be  written  as 

Mvrw  =  2n  ^Wrw  —  i5RwH£la7rZ1  —  fcflA.v|/3 
-  kLx%J  cos  al3  -  kDx sin  al3  +  kLy™pI2\ 

+  ('Wrw  +  °rw )2[£z,Az|Z4  +  kDAxBI6  +  kLx™  cos  al6 
+  kDx$JsmaI6  + kLy%Jl5]}  (47) 


Substituting  the  results  for  the  definite  integrals  I{  — Z6  and 
simplifying  produces 


jyi-B  _  ZAw 
.VRW  —  9 


JcrViUlcosa^Rw  +  ctrw) 


AzB 

4 — (2<Wrw  +  °rW  ~  Vw) 


+  kDJ  i(l)a)RW((^Rw  +  (rRW)[r™  sin  a  +  A.r|] } 


(48) 


The  result  for  the  z  body-axis  is  important  because  it  shows  that 
nonzero  cycle-averaged  forces  can  be  generated  and  used  to  induce 
fore  and  aft  linear  accelerations;  motion  that  would  not  have  been 
possible  using  temporally  symmetric  wingbeat  motion  alone. 

B.  Cycle-Averaged  Moments 

In  this  section,  expressions  for  the  cycle-averaged  aerodynamic 
moments  are  derived.  The  cycle-averaged  moments  are  later  used  in 
the  development  of  a  control  law  that  requires  an  evaluation  of  the 
change  in  cycle-averaged  moments  with  respect  to  a  change  in  the 
control  input  variables. 


Following  a  similar  procedure  for  the  left  wing,  it  can  be  shown  that 

k-LtP LW 


M  B  : 

.VLW 


4pi,(l)  cosa(i5Lw  +  Olw) 


+  — ^  (2ft>LW  +  ctlw  -  <W) 


kDJ ,  (I  )<wLW 


{(Aw  +  ^lw)Kp  sin  a  +  Avf]} 


(49) 


Without  split-cycle  frequency  modulation  there  still  exists  a  nonzero 
cycle-averaged  pitching  moment  if  Az|  =£  Oor  Azf  ^  0.  This  result 
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yields  further  insight  into  vehicle  design  for  controllability.  It 
suggests  that  the  wing  root  hinge  point  should  be  placed  such  that  its 
z  body  location  is  coincident  with  that  of  the  nominal  vehicle  center 
of  gravity,  i.e.,  A z|  =  0,  Az£  =  0.  Such  an  arrangement  will  yield 
a  zero  cycle-averaged  pitching  moment  when  the  split-cycle 
parameters  are  zero,  which  is  a  desirable  feature  for  maintaining 
hover.  Note  that  the  current  Harvard  RoboFly  aircraft  does  not 
incorporate  this  feature  and  would  immediately  pitch  forward  if  not 
constrained  by  the  guide  wires  used  in  their  experiment  [13].  Later,  a 
bob-weight  actuator  is  introduced  that  can  be  used  to  manipulate 
pitching  moments  using  mechanical,  rather  than  aerodynamic, 
means. 


3.  Yawing  Moment 

Expressions  for  the  cycle-averaged  aerodynamic  moments  in  the 
z  body  frame  will  now  be  derived.  Substituting  the  expression  for  the 
instantaneous  z  body  moment  into  Eq.  (10)  yields 


TT-  H  WRW 


A#-Sw  = 


27 r 


jf 


"RW-'SRW 


^zRwuW  dt  + 


/: 


M. 


:RWD 


(t)dt 


“'RW-^RW 


(50) 


Substituting  Eqs.  (3),  (4),  (7),  (8),  (16),  (21),  and  (22)  into  Eq.  (50) 
and  noting  that  the  integrals  are  of  the  form  given  in  Eqs.  (25-30), 
Eq.  (50)  can  be  written  as 


MB  = 


2tt 


( 


(^rw  —  ^rw)“  ~kL  2  1 1  —  kDAxBI2  —  kLy^ppI2 


—  kLx^  cos  al2  —  kDx^y  sinct/i 


+  (wrw  +  °rw)" 


t  wi 
~h2  4 


+  kDAxp5-kLy'^’I6  +  kLx^’cosoiI5  +  kDx^smaI5  |  (51) 


Substituting  the  results  for  the  definite  integrals  I{  — 16  and 
simplifying  yields 


RW  " 


■  i  ( 1)  +  — ^  (2  &>RW  —  5rw  +  (TrW) 


(52) 


geometry,  lift  coefficient,  and  mass  have  been  published  [9,12]  and 
the  expression  derived  by  Doman  et  al.  [12]  is  used  here: 


\jpIAC,Xa)  (54) 

This  condition  assumes  that  there  are  two  wings  flapping 
symmetrically  left  to  right,  i.e.,  Wrw  =  <MlW  and  SRW  =  6LW  =  0, 
and  that  the  angular  position  of  each  wing  in  the  stroke  plane  is 
defined  by  tp(t)  =  cos {coat).  Also,  it  is  assumed  that  the  nominal 
center  of  gravity  of  the  vehicle  and  wing  root  hinges  are  aligned 
such  that  A z£  =  Az£  =  0.  The  control  derivatives  about  hover 
are  presented  in  Table  3.  Note  that  for  Azf  =  Azf  =  0, 

d^.vRW  I  _  n  j,nri  d^fvLW  I  _  n 

dftjRW  ‘fyRW=Wo^RW=0  u  dllU  3 l"LW=(yo,3LW=0  U‘ 


A.  Bob-Weight  Control  Derivatives 

Using  only  the  aerodynamic  control  input  parameters,  changes  in 
pitching  moment  are  accompanied  by  undesirable  changes  in  other 
forces  and  moments.  It  is  for  this  reason  that  a  bob-weight  actuator 
has  been  introduced.  The  bob-weight  should  lie  in  the  vehicle  plane 
of  symmetry  and  at  zero  deflection  should  be  coincident  with  the 
undeflected  center  of  gravity.  The  bob-weight  need  not  be  dead 
weight.  Any  necessary  vehicle  component  that  is  not  sensitive  to 
movement  could  serve  the  purpose,  e.g.,  battery,  power  converter, 
etc.  The  weight  would  be  mounted  near  the  tip  of  a  bimorph 
piezoelectric  actuator.  The  tip  deflection  would  primarily  consist 
of  translation  in  the  z  body-axis  direction,  while  translation  in  the 
x  body  direction  would  be  so  small  as  to  be  negligible.  Hence,  the 
bob-weight  motion  is  in  the  xB,  zfi  plane. 

Let  ma  be  the  total  vehicle  mass  exclusive  of  the  bob-weight  and 
bob-weight  actuator  and  xPgo ,  0,  z£go  be  the  location  of  the  center  of 
mass,  in  body  coordinates,  of  only  the  vehicle.  Let  mBW  denote  the 
mass  of  the  bob-weight  and  actuator  assembly  and  -VBW’  0,  zfw  be  the 
location  of  the  center  of  mass  of  this  assembly  in  body  coordinates. 
Because  of  the  nature  of  the  tip  motion  of  the  bimorph  actuator,  it  is 
assumed  that  motion  of  the  bob-weight  does  not  alter  the  x  location 
of  the  center  of  mass  of  the  bob-weight  assembly;  hence,  vBW  is  a 


Following  a  similar  procedure  for  the  left  wing,  it  can  be  shown  that 


Table  3  Aerodynamic  control  derivatives  about  hover 


JfB 

mz  LW  —  0 


+  ^(2o)RW  —  (SRW  +  crRW) 


(53) 


Without  split-cycle  frequency  modulation,  there  exists  a  nonzero 
cycle-averaged  yawing  moment  on  each  wing  since  y^p  ^  0  and 
uj  0.  When  the  split-cycle  parameters  for  each  wing  are  zero,  the 
cycle-averaged  moments  are  opposing  and  balance  one  another. 
Since  the  fundamental  or  symmetric  wingbeat  frequency  can  be 
independently  varied  for  each  wing,  yawing  moments  can  be 
generated  without  varying  the  split-cycle  parameters. 


Control  derivative  |*=„sf4=o 


1 

_  a^Lw 

1 

95rw 

1  (Urw  —<l)0 

>^rw— 0  95lw 

1  coLW  —o)0  ,dLW  —0 

d^vRW 

1 

_  d^tLW 

1 

9<urw 

1  (Urw  —0)0 

>^rw  =0  9o)Lw 

1  CUlW  —Mo  ’"LW 

d^j’RW 

1 

=  Wyl w 

1 

9<uRW 

1  (Urw  =(O0 

-9rw  — 0  9u>lw 

1  wLW  =a)o  >5lW  =0 

=  0 

=  kLa>„ 
=  0 


M 

loiRW=£U, 


_  | 

IV  =°  95, w  la. 


i  =  o 


(I  ){>>a 


| 

SS,ai  I 


MLW=(WO.Aw=0 


~kDJ,(  l)m0 


V.  Control  Derivatives 

The  parameters  used  to  control  the  aerodynamic  force  and  moment 
profiles  are  the  fundamental  wingbeat  frequencies  wRW  and  <uLW 
and  the  split-cycle  parameters  <5RW  and  i5lw.  For  controllability 
analysis  and  control  synthesis,  the  sensitivity  of  each  cycle-averaged 
generalized  force  to  each  control  input  parameter  is  determined. 

The  partial  derivatives  of  each  cycle-averaged  force  and  moment 
with  respect  to  each  aerodynamic  control  input  parameter  can  be 
calculated  by  differentiating  Eqs.  (33),  (34),  (36),  (37),  (40),  (41), 
(44),  (45),  (48),  (49),  (52),  and  (53)  and  the  full  expressions  for 
these  control  derivatives  are  presented  in  Oppenheimer  et  al.  [19]. 
Here,  the  control  of  the  proposed  aircraft  in  the  vicinity  of  hover  is 
of  interest;  therefore,  the  only  the  control  derivatives  at  the 
hover  condition  where  <wRW  =  ft>LW  =  a>c  and  <5RW  =  <5LW  =  0  are 
presented.  Expressions  relating  the  hover  frequency  to  the  wing 


=<»„.«Ew=o  =  -\kDUo(y%P  +  wT,(l)) 

_  d^rLW  I  _  n 

="o^RW=0  9o»lW  Icvlw=<Wo  Aw=0 

=<d0,jlw=o  =  UD®0(y^p  +  w-Mi)) 

=a>„,«RW=o  =  71(l)ft)0[kixWpcosQ!  +  kD(xf/ 

=(b„,sew=o  =  kLAzBa>0 

=a>„A.w=o  =  JiiX)0)o[kLx^  cost*  +  kD(x™? 

=»)o.4w=o  =  kLAzBLw„ 

—  d^z£w  |  _  n 

=o)0  ,5rw  =0  95lw  I  6>lw  =0)o  >^lw =0 

=r)„,«rw=o  =  — 2<:ift;0[y^p/1(l)  +  “] 
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constant.  Defining  xfg,  0,  zBg  as  the  location  of  the  center  of  mass  of 
the  entire  vehicle,  including  the  bob-weight,  this  assumption  implies 
that  xfg  is  fixed  for  a  given  bob- weight  mass.  Wood  [13]  indicates  that 
the  mass  of  each  wing  is  less  than  1%  of  the  gross  vehicle  weight; 
therefore,  it  is  also  assumed  that  the  mass  of  each  wing  is  negligible. 
Under  the  above  conditions,  the  x  and  z  body  locations  of  the  center 
of  mass  are  given  by 


It  is  assumed  that  a  linear  relationship  exists  between  the  bob-weight 
control  voltage  applied  to  the  actuator,  VbW,  and  the  bob-weight  tip 
deflection: 


SVbw 


—  Aw 


(62) 


x 


B  _ 

eg 


4  om0  +  xgwmBW 
mD  +  mBW 


Additionally,  xfw  is  taken  to  be  constant;  therefore,  =  0;  thus, 
the  control  effectiveness  of  the  bob-weight  can  be  expressed  as 


Zcg  „mo  +  CBWmBW 

ma  +  mBW 


(56) 


9MB  _  9McBg  9zgw 
9VbW  9zbw  9^bw 


(63) 


It  is  important  to  note  that  the  moments  produced  by  aerodynamic 
forces  are  taken  about  the  body  frame  and  are  independent  of  the 
center-of-mass  location,  even  when  the  center  of  mass  is  not 
coincident  with  the  origin  of  the  body-axes  system.  Therefore,  the 
moments  about  the  body  axes  are  composed  of  aerodynamic 
moments  and  moments  due  to  center-of-mass  offset.  In  other  words, 
the  total  cycle-averaged  moments  about  the  body-axes  coordinate 
system  MB  become 

Mb  =  mL„+M*  (57) 

where  MBero  is  the  cycle-averaged  moment  produced  by  lift  and  drag 
of  the  wings  and  M®  is  the  cycle-averaged  moment  due  to  center-of- 
mass  offset  from  the  origin  of  the  body  axes.  The  cycle-averaged 
moment  due  to  center-of-mass  offset  is 

<4  /"rfgXF"  d?  (58) 

■27T  Jo 

where  the  integral  of  each  component  of  the  integrand  is  evaluated, 
rcg  =  [  xcg  0  Zcg  ]  is  the  position  vector  from  the  center  of  mass  to 
the  origin  of  the  body  axes,  and  xfg,  zBg  are  provided  in  Eqs.  (55)  and 
(56).  The  force  in  Eq.  (58)  is  gravitational  and  can  be  written  as 

FcBg  =  RB[-(m0  +  mBW)<?  0  °1  (59) 

where  Rf  is  a  rotation  matrix  from  the  inertial  frame  to  the  body 
frame  and  g  is  the  acceleration  due  to  gravity.  Let  this  force,  in 
general,  be  represented  as  Ffg  =  [Ffcg  Ffcg  Ffcg  ].  Then,  by 


Substituting  Eq.  (56)  into  Eq.  (61)  and  computing  the  derivatives  as 
specified  in  Eq.  (63)  yields 


9M“ 

9  Try 


-Flk* 


m0+n ibw  *  BW 

”BW  Ff  £bw 

mo  i  mBW  xeg  DVV 


(64) 


For  flight  near  hover,  Ff  Fy  because  both  the  vehicle  pitch  and 
yaw  attitudes  stay  near  0  deg. 


B.  Control  Effectiveness 

The  final  form  of  the  control-oriented  model  consists  of  a  cycle- 
averaged  blade-element  model  that  is  linearized  about  hover.  Within 
this  framework,  a  linear  relationship  exists  between  the  cycle- 
averaged  force  and  moment  commands  generated  by  a  conventional 
feedback  control  law  and  the  control  parameters.  Letting  G  denote  a 
generalized  force  or  moment,  the  total  increment  of  commanded  G  is 


99, 


9G  „  9G  A  9G  9G  A 

<JRW  +  q - 4%w  +  TCg - Aw  +  - - Awlw 


9ft>R 


99t 


9w, 


9G 

+  wr, —  Tbw 

9Krw 


(65) 


Writing  this  expression  for  the  three  forces  and  three  moments 
produces 


'  A  Fjom 
AFr 

_ )fcom 

A  F  B 

z-i_Vzcom 

A  M  B 

L-*1Y1xcom 

A  WB 

_ ycom 

AAFS ... 


0  Bi  2  0 
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~kDJ  i  (I  )0->o  0  —kDJi(l)co„ 

~\kD0)o(y^  +  wJ  i(  1))  0  \kDa>0(y™  +  wJ^l)) 

B5i  0  B53 

0  b62  0 


6,4 

0 

0 

0 

0 

B(,4 


0 

0 

0 

b45 

b55 

0 


Aw 

Aw 

Awrw 

A  wRW 

'Aw 

=  BA 

Aw 

A®lw 

Awlw 

^BW 

'/BW 

(66) 


performing  the  cross-product  operation,  the  moment  due  to  center- 
of-mass  offset  is 


B  pB 

^Cg*  2C-~ 


-7B  FB 

ycg 


■  4  FB 

Ccg 


rB 

'lcgJ  rcg 


df 


(60) 


where  BA  is  the  control  effectiveness  matrix.  The  elements  of  this 
matrix  are  the  sensitivities  of  the  cycle-averaged  forces  and  moments 
to  changes  in  control  input  variations  that  are  presented  in  Table  3  and 
Eq.  (64).  Also, 


B\2—B\n  —  kLu>0  B62  —  —BM  —  —2kLa>c 


4P-/  id)+- 


B5I  =  7i  ( 1  )(D0[kLxT’  cosot  +  kD(x™p  sina  +  Axf )] 


If  the  position  vector  rBg  and  the  force  in  Eq.  (59)  are  constant  over  a 
full  wingbeat  cycle,  then  Eq.  (60)  becomes 


M 


B  _ 
eg  — 


~,B 


FB 


ycg 

7B  fB  _  XB  r?B 

t-cg1  Xcg  -^cg1  Zcg 

xB  Fb 

eg  ycg 


(61) 


B53  =  J  i  ( 1  )co0[kLx™  cos  a  +  kD(x*j?  sina  +  Axf )] 


Bic  — 


,  (  mBW  A 

—  ~kBW  I - ■ - I 

\m0  +  mBWJ 


FB 

rcg 


cp 

B55  =  A 


(  '”BW  \fb 

\m0  +  mBWJ 


(67) 


A  feedback  control  law,  discussed  shortly,  generates  the  cycle- 
averaged  force  and  moment  commands  that  form  the  left-hand  side  of 
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Eq.  (66).  The  control  allocation  objective  is  to  find  the  control  input 
vector  [  i5rw  A&>rw  <5lw  A&>lW  Vbw  ]r,  such  that  Eq.  (66)  is 
satisfied. 

Examination  of  this  matrix  reveals  what  can  be  achieved  by 
manipulating  the  four  parameters  that  define  the  motion  of  the  wings 
and  the  bob-weight.  The  rank  of  this  matrix  is  five.  Most  obviously,  it 
is  not  possible  to  generate  direct  side  forces  on  the  vehicle  using  the 
five  control  parameters.  The  next  observation  is  that  rows  associated 
with  the  z  force  and  pitching  moment  would  be  linearly  dependent 
for  Az|  =  A;f  =  OandA*!  =  Avf  if  it  were  not  for  the  presence  of 
the  bob-weight  actuator.  These  conditions  would  exist  on  a  well 
designed  vehicle  because  of  the  existence  of  a  plane  of  symmetry  and 
the  fact  that,  from  a  vehicle  design  perspective,  it  is  desirable  that  the 
z  offset  between  the  nominal  center  of  gravity  and  the  wing  root 
hinge  points  be  zero.  Such  an  arrangement  eliminates  a  source  of 
undesirable  pitching  moments  associated  with  a>0  that  would  tend  to 
destabilize  the  hover  condition  unless  active  control  were  applied. 
The  sparsity  of  the  control  effectiveness  matrix  shows  that  the  split- 
cycle  method  provides  a  high  degree  of  input  decoupling,  which 
simplifies  control  allocation.  The  analytic  formulation  allows  for  a 
great  deal  of  insight  into  the  physics  of  the  control  problem,  which 
can  aid  in  the  design  of  the  controller  as  well  as  the  vehicle  itself  [19]. 


VI.  Equations  of  Motion 

A  standard  set  of  6-degree-of-freedom  rigid  body  equations  of 
motion  are  employed  [22],  These  are  the  same  equations  that 
describe  the  motion  of  a  rigid  body  aircraft;  however,  because  the 
angle  of  each  wing  in  the  stroke  plane  is  a  periodic  function  of  time, 
the  forces  and  moments  that  drive  the  equations  consist  of  multiple, 
independent,  and  variable  periodic  functions  of  time.  It  is  assumed 
that  the  wings  are  massless;  therefore,  there  exists  no  inertial 
coupling  between  the  fuselage  and  the  wings.  It  is  hypothesized  that 
the  inertial  coupling  between  the  wings  and  body  will  be  small 
because  in  the  Harvard  Robofly  the  wings  constitute  less  than  1%  of 
the  total  vehicle  weight  [13].  The  equations  of  motion  are 
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where  I  is  the  inertia  matrix;  to  =  [p  q  r]  is  the  angular  rate 
vector;  M B(t)  is  the  sum  of  the  left  and  right  wing  instantaneous 
aerodynamic  moment  vectors  consisting  of  the  rolling,  pitching,  and 
yawing  moments;  u,  v ,  and  w  are  the  translational  velocities;  F B(t)  is 
sum  of  the  left  and  right  wing  instantaneous  aerodynamic  force 
vectors;  m  is  the  vehicle  mass;  g  is  the  acceleration  due  to  gravity;  x, 
y,  and  z  are  the  vehicle  positions  with  respect  to  an  inertial  frame;  and 
q  =  <jr0  +  q{i  +  q2j  +  q2k  is  a  quaternion  used  to  perform  the  3-2-1 
standard  Euler  transformation.  Also,  Rg  =  Rf  are  rotation  matrices 
that  transform  from  body  to  inertial  axes  and  back. 


A.  Blade-Element  Simulation 

A  low-fidelity  blade-element-based  model  of  the  instantaneous 
aerodynamic  forces  and  moments  is  used  to  drive  equations  of 


motion  in  the  first  simulation  model.  These  instantaneous  forces  and 
moments  were  presented  earlier  in  Table  1  and  Eqs.  (21-24).  This 
model  does  not  account  for  leading-edge  vortex  and  wake-capture 
effects  that  have  been  observed  in  experiments.  The  purpose  of  this 
simulation  model  is  to  assess  the  efficacy  of  the  cycle-averaged 
control  law  when  applied  to  a  vehicle  model  that  includes  time 
periodic  variations  in  the  forces  and  moments  due  to  the  flapping 
wings. 


B.  Unsteady  Aerodynamic  Simulation 

Unsteady  aerodynamic  phenomena  can  have  a  significant 
influence  upon  the  instantaneous  forces  and  moments  produced  by 
flapping  wings  at  low  Reynolds  numbers.  Significant  discrepancies 
between  the  instantaneous  aerodynamic  forces  predicted  by  blade- 
element  theory  and  experimental  measurements  have  been  observed; 
however,  it  has  been  found  that  in  a  cycle-averaged  sense,  blade- 
element  theory  provides  a  reasonable  prediction  of  cycle-averaged 
lift  [20].  Blade-element  theory  has  been  found  to  underpredict  cycle- 
averaged  drag  when  compared  with  experimental  measurements. 
Since  the  cycle-averaged  control  law  is  based  on  the  sensitivities  of 
cycle-averaged  lift  and  drag  with  respect  to  the  adjustable  wing 
kinematic  parameters  and  not  the  actual  magnitudes  of  the  forces,  it  is 
natural  to  ask  whether  the  cycle-averaged  blade-element  model  used 
to  derive  the  control  laws  captures  these  sensitivities  with  sufficient 
fidelity  to  provide  a  reasonable  level  of  closed-loop  stability  and 
tracking  performance.  To  address  this  question,  the  fidelity  of  the 
aerodynamic  forces  produced  by  the  flapping-wing  vehicle  model  is 
increased  by  making  use  of  the  experimental  results  produced  by 
Sane  and  Dickinson  [20].  The  experimental  results  were  used  to 
produce  time-varying  lift-  and  drag-coefficient  profiles  that  provide  a 
more  complete  representation  of  the  forces  and  moments  produced 
by  a  flapping  wing  at  low  Reynolds  numbers.  This  representation  not 
only  includes  translational  and  rotational  contributions  to  the  forces 
that  can  be  predicted  by  quasi-steady  methods,  but  also  includes 
contributions  from  unsteady  leading-edge  vortex  and  wake  inter¬ 
action  effects. 

The  lift  and  drag  profiles  for  a  given  stroke  and  angle-of-attack 
profile  are  extracted  from  measurement  time  histories.  The  cor¬ 
responding  instantaneous  aerodynamic  coefficients  were  computed 
to  form  a  look-up  table.  The  computation  of  the  coefficient  tables  is 
performed  using  the  parameters  of  the  flapping-wing  system 
employed  by  Sane  and  Dickinson  [20],  while  the  computations  of  the 
aerodynamic  forces  in  the  simulation  model  used  for  this  study  use 
the  parameters  of  the  proposed  flapping-wing  MAV.  The  lift  and  drag 
force  measured  in  the  Sane  and  Dickinson  experiment  were  observed 
for  a  90  deg  stroke  amplitude  with  a  symmetric  flip  of  the  wing,  about 
the  spar  axis,  that  occurred  in  the  vicinity  of  wing  stroke  reversals 
[20].  The  measurement-based  data  was  approximated  by  producing 
1000-point  look-up  tables  for  lift  and  drag  coefficients  and  for  angle 
of  attack,  using  constant  interval  sampling  of  the  phase  of  the  stroke. 
The  phase  of  the  stroke  waveform  is  used  as  an  index  to  extract  data 
from  the  resulting  aerodynamic  coefficient  tables.  The  angle  of  attack 
is  also  obtained  from  a  look-up  table  that  corresponds  to  the  set  of 
force  coefficients  being  used.  The  angle-of-attack  table  is  necessary 
because  the  instantaneous  values  of  the  aerodynamic  coefficients 
depend  upon  the  instantaneous  angle-of-attack  profile.  The 
appropriate  rotational  angle  is  therefore  produced  to  transform  the 
lift  and  drag  forces  in  the  spar  frame  to  the  body  frame.  A  comparison 
of  the  unsteady  table  look-up  model  to  the  instantaneous  blade- 
element  and  cycle-averaged  blade-element  models  for  a  sinusoidal 
velocity  profile,  i.e.  &>RW  =  &>LW  =  a>0  andi5RW  =  5LW  =  0,  is  shown 
in  Fig.  5.  It  can  be  seen  that  the  instantaneous  blade-element  model 
underpredicts  drag  but  captures  the  magnitude  of  the  instantaneous 
lift  rather  well.  The  cycle-averaged  blade-element  model  used  for  the 
model-based  control  design  differs  significantly  from  both  the 
instantaneous  blade-element  and  unsteady  aerodynamic  model  over 
the  course  of  the  wingbeat;  however,  the  controller  is  designed  to 
control  the  low-frequency  behavior  of  the  fuselage,  which  is  well 
below  the  wingbeat  frequency. 
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Fig.  5  Aerodynamic  model  comparison  for  temporally  symmetric  wingbeat. 


In  summary,  the  unsteady  model  relaxes  the  following 
assumptions  used  in  the  derivation  of  the  control-oriented  model: 
the  two-dimensional  sectional  aerodynamic  coefficients  are  known 
and  constant  throughout  each  stroke,  and  the  passive  wing  rotation 
joint  is  on  a  limit  when  the  wing  angular  velocity  in  the  stroke  plane  is 
nonzero.  The  first  assumption  is  replaced  with  a  model  where  the 
aerodynamic  coefficients  of  each  wing  are  known,  but  vary  as  a 
function  of  angular  position  in  the  stroke  plane.  The  second  assump¬ 
tion  is  replaced  by  a  model  where  the  angle  of  attack  smoothly  varies 
as  the  wing  moves  through  the  stroke  plane. 


VII.  Inner-Loop  Control 

The  objective  of  the  inner-loop  controller  is  to  render  a 
neighborhood  of  a  given  point  and  orientation  attractive.  Seeking 
zero-error  tracking  to  a  fixed  point  and  orientation  is  not  possible 
because  of  the  inherent  periodic  forcing  of  the  system  that  results 
from  the  flapping  wings.  The  control  strategy  is  to  use  the  available 
analytical  expressions  for  wingbeat  cycle-averaged  forces  and 
moments  as  a  control-oriented  model  and  consider  the  averaged 
forces  and  moments  as  virtual  control  inputs.  The  controller, 
therefore,  assigns  the  wingbeat  parameters  for  each  cycle;  thus,  the 
output  of  the  controller  need  only  be  computed  once  per  wingbeat. 
Analysis  of  the  cycle-averaged  forces  and  moments  shows  that  using 
a  combination  of  left  and  right  fundamental  frequencies,  split-cycle 
parameters,  and  bob-weight  shifting  yields  control  authority  over  all 
but  the  cycle-averaged  side  force.  Hence,  the  controller  will  be  a  five- 
DOF  controller,  but  it  is  expected  to  perform  six-DOF  tracking  tasks. 

The  position,  velocity,  attitude,  and  angular  velocities  of  the  body 
are  assumed  to  be  available  for  feedback.  They  are  used  to  compute 
an  error  between  the  vehicle  state  and  a  desired  point  and  orientation. 
The  high  degree  of  input  decoupling  present  in  the  control  matrix  BA 
enables  the  use  of  a  simplistic  feedback  gain  structure  where  the 
approach  is  to  enforce  a  second-order  response  in  each  degree  of 
freedom,  considering  each  force  and  moment  separately.  The  output 
of  the  feedback  gains  is  the  desired  cycle-averaged  forces  and 
moments  and  forms  the  left-hand  side  of  Eq.  (66).  Finally,  control 
allocation  is  used  to  compute  the  wingbeat  parameters,  where  the 
allocation  is  based  on  the  control  derivatives  derived  from  the 


analytic  expressions  for  the  wingbeat  cycle-averaged  forces  and 
moments. 

A.  Error  Computation 

The  error  computation  employs  a  given  reference  position  and 
orientation,  expressed  with  respect  to  the  inertial  coordinate  frame. 
The  error  is  expressed  in  body  coordinates,  transformed  appro¬ 
priately  from  the  difference  between  the  reference  and  the  state  of  the 
vehicle.  The  feedback  gain  structure  is  based  on  considering  each 
degree  of  freedom  separately  in  the  body  coordinate  frame,  because 
expressions  for  the  average  forces  and  moments  are  available  in  the 
body  frame.  Computing  the  error  in  body  frame  coordinates  is 
therefore  consistent  with  the  desired  feedback  gain  structure.  The 
velocity  error  can  be  computed  similarly,  when  the  reference  is 
allowed  to  vary.  The  angular  velocity  error  is  approximated  by  the 
angular  velocity  of  the  vehicle  in  body  coordinates.  For  a  fixed 
reference  orientation,  no  approximation  is  needed. 

A  desired  position  xd  and  a  desired  orientation  defined  by  the  Euler 
angles  <I>f  are  given,  where  D  denotes  the  desired  body  frame. 
Alternatively,  the  reference  orientation  is  represented  by  the  rotation 
matrix  Rf  from  the  inertial  frame  to  the  desired  frame.  The  tracking 
error  is  computed  as 

xf  =  R  f(x'-x'd)  vf=xB-Rf(x'd)  R®  =  Rf  (Rf)r 

<°e=(P  1  rY  (72) 

where  xB  =  (xe  ye  ze)T  is  the  position  error,  vB= 
(if  yB  zB)T  is  the  velocity  error,  R£  is  the  rotation  matrix 
from  the  desired  frame  to  the  body  frame,  and  =  ( pe  qe  re)T 
is  the  angular  velocity  error.  A  vector  of  Euler  angle  tracking  errors 
[4>e,  8e,  ];  can  be  derived  from  the  R£  rotation  matrix  using 

standard  methods  [23]. 

B.  Feedback  Gain  Structure 

The  feedback  controller  is  based  upon  the  simplistic  view  of 
treating  each  degree  of  freedom  separately.  This  is  possible  because  of 
the  use  of  a  control  allocator  that  translates  cycle-averaged  gener¬ 
alized  force  commands  into  wingbeat  oscillator  and  bob-weight 
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Fig.  6  Five-degree-of-freedom  split-cycle  controller. 
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actuator  commands.  Additionally,  near  the  hover  condition,  the 
equations  of  motion  reduce  to  a  bank  of  decoupled  double  integrators, 
because  there  are  no  aerodynamic  forces  acting  on  the  fuselage  other 
than  those  produced  by  the  beating  wings,  which  are  accounted  for  in 
the  control  effectiveness  matrix.  Thus,  each  axis  is  considered  as  a 
simple  second-order  one-dimensional  system: 

X  =  u  (73) 

Defining  Xd  as  a  desired  reference  command,  it  is  well-known  that 
assigning  u  =  —u>2(x  ~  Xd)  ~  2£o>„(x  —  Xd)  forces  the  response  of 
the  system  to  that  of  a  second-order  damped  harmonic  oscillator  with 
natural  frequency  a>n  and  damping  ratio  £ .  This  simple  view  is  applied 
to  each  degree  of  freedom  of  the  system.  A  feedback  structure  is 
assigned  to  the  commanded  cycle-averaged  forces  and  moments  as 
follows: 

i;Fxa)Fxn)m 

^L,  =  (-Z^Fzn  -  ie^Fz(0Fzn)m 
=  (~<Pe0)2Mxn  -  pe2^MxcoMxn)Ixx 

Mycom  (  ^e^Myn  ^1  e2^My^Myn)^yy 

^Zcom  ^  tye^Mzn  ^e2^Mz^Mzn)^zz  (74) 

A  block  diagram  of  the  control  structure  is  shown  in  Fig.  6.  An 
alternative  view  is  to  consider  the  structure  as  a  sequential  loop 


closure  controller  that  uses  proportional-derivative  feedback,  where 
the  gains  are  tuned  by  choosing  appropriate  values  for  each  a>.„  and  f. 
in  Eq.  (74). 

Note  that  the  v-axis  translation  is  not  directly  controlled,  in  the 
sense  that  there  is  no  Ff  .  Instead,  the  y-axis  translation  is  regulated 
by  injecting  ye  into  the  desired  yaw  assigned  by  an  outer-loop 
controller,  which  is  described  in  the  subsequent  section.  Assuming  a 
nominal  Fx  cycle-averaged  force  equal  to  weight,  mg,  and  zero  side 
force,  then  elementary  trigonometry  shows  that  at  zero  pitch  and 
roll,  the  force  directed  along  the  y  inertial  frame  axis  is  sin(i jr)mg,  or 
\[smg  with  a  small  angle  approximation.  This  simplified  analysis  is 
employed  to  construct  a  desired  approximate  side  translation  force, 
through  an  added  yaw  command.  The  added  desired  yaw  has  the 
same  feedback  gain  structure  as  the  inner-loop  feedback,  and  is 
given  by 

f  d  =  (-yBeO)2yn  -  v?2 i;y(0yn)  -  (75) 

§ 

The  y  translation  reference  command  i //d  appears  in  the  outer-loop 
controller  presentation. 

C.  Control  Allocation 

A  pseudo-inverse  control  allocation  scheme  [24]  is  used  to 
determine  the  necessary  control  inputs  to  produce  the  desired  cycle- 
averaged  forces  and  moments.  The  control  allocator  computes  the 
fundamental  frequencies,  split-cycle  parameters,  and  voltage  applied 


Table  4  Properties  of  proposed  flapping- wing  MAV 


Fig.  7  Desired  trajectory. 


Variable 

Value 

Units 

m„ 

60 

mg 

mbw 

20 

mg 

Height 

11 

mm 

Width 

4 

mm 

Depth 

1 

mm 

Ia 

1395 

mm4 

^BW 

0. 165e  —  5 

m 

Volt 

Nominal  CG  location 

[5.5  0  0] 

mm 

^max 

4 

mm 

d„ 

3 

mm 

R 

15 

mm 

Arf 

[3.5  2  0] 

mm 

Arf 

[3.5  -2  0] 

mm 

«£/ 

45 

deg 

aD 

45 

deg 

Hover  frequency 

113.61 

Hz 

=o=l  mg 

y  .5pCL(a)IA 
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Table  5  Inner-loop  and  y-body-axis 
regulation  controller  parameters 


Tuning  parameter 

Value 

<°Fxn 

7 

%Fxn 

1 

& Fyn 

7 

£ Fyn 

1 

toFzn 

4 

Gzn 

.85 

to\fxn 

40 

%Mxn 

1 

toMyn 

98 

£ Myn 

2.7 

tofAzn 

30 

KftAzn 

1.1 

to  the  bob-weight  actuator  based  on  the  cycle-averaged  force  and 
moment  commands.  The  fundamental  frequencies  and  split-cycle 
parameters  govern  the  behavior  of  oscillators  that  generate  time- 
varying  wing  position  commands.  Inputs  to  the  oscillators  are  held 
constant  over  the  course  of  a  wingbeat  in  order  to  avoid  changing  the 
shape  of  the  split-cycle  waveform  in  midcycle.  The  element  used  to 
accomplish  this  function  is  called  a  cycle  zero-order  hold.  The  output 
of  the  oscillators  are  applied  to  the  bimorph  piezoelectric  actuators.  It 
is  assumed  that  the  actuator  bandwidth  is  sufficient  to  track  the  wing 
position  commands. 

Given  the  control  variable  vector  u= 
[i5rw  Awrw  5lw  Aa>LW  VBW  ]7  and  hover  condition  ua  = 
[0.  a>0 ,  0,  a>0,  0]7,  the  control  allocation  at  time  tk  is 

expressed  as  follows: 

«  (4)  =  B\Gcom{tk )  +  ua  (76) 


operations  is  relaxed  by  replacing  the  demand  of  perfect  alignment  of 
the  z  body  axis  with  requiring  that  the  z  axis  be  within  (ft  of  being 
directed  towards  the  waypoint,  creating  a  2(/>-wide  heading  window. 
Furthermore,  once  the  projection  of  the  vehicle’s  position  on  the 
(y,  z)  plane  of  the  inertial  frame  is  within  distance  x  of  the  waypoint, 
the  desired  orientation  is  not  allowed  to  change,  creating  a  position 
window  of  radius  x.  The  position  window  prevents  the  outer-loop 
controller  from  perpetually  prescribing  oscillating  roll  commands  as 
the  vehicle  cycles  about  the  desired  waypoint. 

Without  loss  of  generality,  it  is  assumed  that  the  initial  position 
before  a  waypoint  change  is  zero  and  that  the  vehicle  is  hovering; 
i.e.,  the  x  axis  of  the  body  frame  is  vertical.  Given  waypoint 
=  ( xwp  vwp  zwp  )r,  the  waypoint  orientation  is  computed  as 


(arctan(— (ywp  -  y),  (zwp  -  z))  \ 

o  j  ™ 

wherex7  =  [x  y  z  ]r  is  the  position  of  the  vehicle  and  </>wp  is  the 
desired  heading  derived  from  elementary  trigonometry,  expressed  in 
terms  of  a  four-quadrant  inverse  tangent.  The  vector  denotes  the 
additional  attitude  command  generated  by  feedback.  This  vector 
consists  only  of  i j/d  from  Eq.  (75),  but  additional  commands,  such  as 
a  desired  pitch,  can  be  generated  to  assist  forward/backward 
translation. 

The  desired  orientation  is  generated  as  follows: 


G„(*?)  +  <*>? 
G„(*?)  +  <*>? 


if  ^Cy  -  TwP)2  +  (Z  -  4vp)2  5  X 

otherwise 


(79) 


where 

Gcom  =  [Flm  0  F?com  MBXam  Myc<m  M?com]T  (77) 

and  BA  denotes  the  pseudo-inverse  of  BA.  The  structure  of  the 
control  allocation  matrix  is  simplified  because  the  split-cycle  control 
strategy  yields  a  highly  decoupled  control  effectiveness  matrix.  Note 
that  the  term  FB  is  identically  set  to  zero,  since  direct  side  force 
cannot  be  generated  using  the  split-cycle  method  and  that  side-force 
corrections  are  made  by  adjustments  to  yaw  attitude. 

VIII.  Outer-Loop  Controller 

The  outer-loop  controller  generates  the  desired  position  and 
orientation,  with  respect  to  the  inertial  coordinate  frame,  which  is  fed 
to  the  error  generator.  This  includes  the  tasks  of  providing  desired 
trajectories  for  waypoint  tracking  and  regulating  the  y-axis 
translation  by  generating  a  desired  yaw,  since  the  inner-loop  does  not 
directly  regulate  y  translation.  The  desired  position  and  orientation 
vary  as  step  changes  are  fed  through  second-order  continuous-time 
smoothing  filters. 


where  Go(0  represents  a  second-order  filter  bank  with  natural 
frequency  and  damping  ratio  tu#c,  respectively,  and  is  the 
value  of  held  from  the  instant  when  the  condition 

y  —  ywp)2  +  (z  —  zwp)2  5  x  was  most  recently  satisfied.  The 

filter  G„  is  used  to  maintain  differential  continuity  of  the  desired 
orientation  after  a  discrete  waypoint  change.  Note  that  the  position 
condition  in  (79)  is  a  condition  on  the  difference  between  the  position 
of  the  vehicle  and  the  waypoint,  projected  on  the  (y,  z)  plane  of  the 
inertial  frame;  it  is  not  a  function  of  the  desired  position  described 
below. 

Because  of  the  body  frame  orientation  at  hover,  it  is  convenient  to 
consider  as  a  1-2-3  rotation.  This  is  comparable  to  the  common  3- 
2- 1  rotation  for  a  standard  aircraft  in  level  flight,  where  the  nose  is 
directed  forward  along  the  x  axis  and  the  z  axis  is  directed  downward, 
while  the  flapping-wing  vehicle  at  hover  is  oriented  such  that  its 
ventral  side  is  directed  forward  along  the  z  axis  and  the  x  axis  is 
directed  upward.  The  construction  of  3>f  and  $7P  is  thus  made 
simple,  with  (ftd  as  the  desired  heading  of  the  vehicle.  The  rotation 
matrix  for  =  [  (ftd,  0d,  \ftd  ]r  is 


Rf  = 


cos  (fd)  cos  (6d)  sin^)  cos(cpd)  +  cos(i/rrf)  sin(6d)  sin  (<ftd)  sm(fd)  sin  (<ftd)  -  cos(^)  sin(£rf)  cos  (<ftd)  \ 

-  sin(i ftd)  cos(dd)  cos(Vg)  cos (<ftd)  -  sinfi/^)  sin(6»rf)  sin(c^)  cos(fd)  sin (cftd)  +  sm(iftd)  sin(£rf)  cos {<pd) 
sin(6>rf)  -  cos(erf)  sin(0rf)  cosf^jcos (<pd)  ) 


The  waypoint  reference  generator  computes  a  desired  trajectory 
for  the  inner-loop  controller  based  on  a  roll-and-translate  approach, 
which  is  similar  to  a  bank-to-tum  strategy.  The  approach  is  to  use  the 
roll  angle  to  first  adjust  the  heading  angle  of  the  hovering  vehicle, 
such  that  the  positive  z  axis  of  the  body  frame  is  directed  towards  the 
waypoint.  Then,  the  aircraft  translates  to  the  waypoint  while  the 
attitude  is  regulated.  Hence,  the  desired  position  is  the  waypoint  itself 
and  the  desired  orientation  has  zero  pitch  and  yaw  attitude,  but  the 
roll  is  determined  by  the  desired  heading.  The  separation  of  these  two 


which  is  used  in  the  error  generation,  Eq.  (72),  to  compute  the  attitude 
error. 

The  desired  position  is  computed  as 


/  _  {  Gp(: 
d~  Up 


(Xwp)  if  |  (<Pd  (ft d)  0wp  A  (ft 

otherwise 


(80) 
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where  Gp(-)  represents  a  second-order  filter  with  natural  frequency 
and  damping  ratio  coxim,  respectively,  and  x^p  is  the  value  of 
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Table  6  Outer-loop  controller 
and  filter  parameters 


Tuning  parameter 

Value 

Gp  natural  frequency 

4 

Gp  damping  ratio 

1 

G0  natural  frequency 

20 

G0  damping  ratio 

1 

X 

5  ■  10~2  m 

4> 

75° 

Gp(x[ ,p),  held  from  the  instance  when  the  condition  \((j>d  —  (f>d)  — 
$wp|  5  0  was  most  recently  satisfied.  The  filter  Gp  is  used  to 
maintain  differential  continuity  of  the  desired  orientation  after  a 
discrete  waypoint  change. 

The  conditions  in  Eqs.  (79)  and  (80),  representing  the  heading  and 
position  windows,  depends  on  the  position  of  the  vehicle,  projected 
on  the  (y,  z)  plane  of  the  inertial  frame.  In  the  case  of  Eq.  (80),  the 
dependency  originates  from  the  waypoint  orientation  Eq.  (78),  which 
depends  on  the  position  of  the  vehicle.  The  window  is  set  such  that 
the  controller  can  be  tuned  to  render  a  subset  of  the  interior  of  the 
window  attractive. 


Fig.  8  Inertial  position  time  history  resulting  from  waypoint  tracking  test  using  blade-element  and  unsteady  simulations. 
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Fig.  9  Attitude  time  history  resulting  from  waypoint  tracking  test  using  blade-element  and  unsteady  simulations. 
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Fig.  10  Time  history  of  fundamental  and  split-cycle  frequencies  resulting  from  waypoint  tracking  test  using  blade-element  and  unsteady  simulations. 


IX.  Results 

Figure  7  shows  a  three-dimensional  trajectory  that  the  vehicle  is 
commanded  to  follow.  The  vehicle  begins  at  the  origin  in  inertial 
space;  waypoints  1  through  5  are  then  fed  to  the  outer-loop  controller 
at  2-s  intervals  in  ascending  order.  Vehicle  properties  are  shown  in 
Table  4.  The  inner-  and  outer-loop  control  parameters  used  to 
generate  the  subsequent  simulation  results  are  provided  in  Tables  5 
and  6,  respectively. 

Figures  8  and  9  show  the  commanded  and  actual  positions  and 
attitudes  of  the  vehicle  center  of  gravity  used  to  achieve  waypoint 
following  resulting  from  the  instantaneous  blade-element  and 


unsteady  aerodynamic  simulations.  In  both  cases,  tracking  per¬ 
formance  is  acceptable  and  the  vehicle  is  capable  of  performing  the 
desired  maneuvers;  however,  the  performance  in  the  unsteady 
aerodynamic  simulation  is  noticeably  degraded  due  to  larger 
modeling  errors.  Figure  10  shows  the  resulting  left  and  right  wing 
fundamental  and  split-cycle  frequencies.  It  can  be  observed  that  the 
left  and  right  wing  fundamental  and  split-cycle  frequencies  are 
significantly  more  dynamic  in  the  unsteady  aerodynamic  simulation 
than  those  that  were  generated  in  the  blade-element  simulation. 
In  both  cases,  differences  between  the  left  and  right  split-cycle 
frequencies  result  in  the  production  of  rolling  and  yawing  moments 


x  itT* 


0  2  4  6  8  10 

Time  (s) 

Fig.  11  Bob-weight  position  and  eg  shift  in  z  body  direction  due  to  bob-weight  movement;  waypoint  tracking  test  using  blade-element  and  unsteady 
simulations. 
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that  are  used  to  align  the  vehicle’s  heading  with  the  next  waypoint 
and  to  regulate  lateral  translation.  The  position  of  the  bob- weight  and 
the  location  of  the  vehicle  center-of-gravity  in  the  z  body-axis 
direction  are  shown  in  Fig.  1 1 .  In  the  blade-element  simulation,  it  can 
be  seen  that  the  bob-weight  shifts  the  center  of  gravity  by  about 
—0.08  mm  in  the  worst  case.  This  corresponds  to  a  —0.32  mm  shift 
in  the  z  position  of  the  bob-weight.  In  the  unsteady  simulation,  it  can 
be  observed  that  bob-weight  activity  is  much  higher.  The  bob-weight 
shifts  the  center  of  gravity  by  up  to  0.6  mm  in  the  worst  case,  which 
corresponds  to  a  2.2  mm  shift  in  the  z  position  of  the  bob-weight.  The 
magnitude  of  the  required  bob-weight  in  the  unsteady  aerodynamic 
simulation  movement  is  considered  excessive  and  might  necessitate 
the  use  of  an  alternative  actuator  because  of  the  limited  range  of 
motion  of  bimorph  piezoelectric  actuators.  Mechanical  amplification 
of  the  bimorph  tip  motion  or  an  alternative  pitch  control  strategy  may 
have  to  be  employed  in  practice. 

X.  Conclusions 

Modeling  and  control  of  a  flapping-wing  MAV  were  considered  in 
this  manuscript.  The  results  support  the  conclusion  that  it  is  possible 
to  independently  control  vertical  and  horizontal  forces  as  well  as 
rolling  and  yawing  moments  using  two  physical  actuators  that  allow 
the  angular  velocity  profiles  of  the  left  and  right  wings  to  vary  from 
upstroke  to  downstroke.  The  addition  of  an  actuated  bob-weight 
enabled  control  of  the  aircraft  pitching  moment.  A  method  called 
split-cycle  constant-period  frequency  modulation  generated  the 
required  wingbeat  actuator  position  commands  and  resulted  in  a  high 
level  of  input  decoupling,  which  enabled  controllers  with  a  simple 
structure  to  successfully  provide  controlled  maneuverability  in  the 
vicinity  of  hover.  The  method  has  the  desirable  feature  that  cycle- 
averaged  forces  and  moments  can  be  analytically  related  to  the 
wingbeat  kinematics  and  vehicle  design  parameters,  a  feature  that 
facilitates  multidisciplinary  design  efforts  where  controllability  must 
be  considered  from  the  outset.  Additionally,  simulation  results  reveal 
that  a  model-based  control  law,  derived  from  a  cycle-averaged  blade- 
element  aerodynamic  model,  can  successfully  control  simulation 
models  that  include  unsteady  aerodynamic  effects. 

References 

[1]  Deng,  X.,  Schenato,  L.,  Wu,  W.  C.,  and  Sastry,  S.  S.,  “Flapping  Flight 
for  Biomimetic  Robotic  Insects,  Part  1 :  System  Modeling,”  IEEE 
Transactions  on  Robotics  and  Automation,  Vol.  22,  No.  4,  2006, 
pp.  776-788. 

[2]  Deng,  X.,  Schenato,  L.,  Wu,  W.  C.,  and  Sastry,  S.  S.,  “Flapping  Flight 
for  Biomimetic  Robotic  Insects,  Part  2:  Flight  Control  Design,”  IEEE 
Transactions  on  Robotics  and  Automation,  Vol.  22,  No.  4,  2006, 
pp.  789-803. 

[3]  Khan,  Z.  A.,  and  Agrawal,  S.  K.,  “Control  of  Longitudinal  Flight 
Dynamics  of  a  Flapping- Wing  Micro  Air  Vehicle  Using  Time-Averaged 
Model  and  Differential  Flatness  Based  Controller.”  Proceedings 
of  the  American  Control  Conference,  IEEE,  New  York,  July  2007, 
pp.  5284-5289. 

[4]  Dileo,  C.,  and  Deng,  X.,  “Design  of  and  Experiments  on  a  Dragonfly- 
Inspired  Robot,”  Advanced  Robotics,  Vol.  23,  Nos.  7-8,  2009, 
pp.  1003-1021. 

doi:  10. 1 163/156855309X443160 

[5]  de  Croon,  G.  C.  FI.  FI.,  de  Clercq,  K.  M.  E.,  Ruijsink,  R.,  Remes,  B.,  and 
de  Wagter,  C.,  “Design,  Aerodynamics,  and  Vision-Based  Control  of 
the  DelFly,”  International  Journal  of  Micro  Air  Vehicles,  Vol.  1 ,  No.  2, 
June  2009,  pp.  71-97. 


[6]  Jones,  K.  D.,  Bradshaw,  C.  J.,  Papadopoulos,  J.,  and  Platzer,  M.  F., 
“Bio-Inspired  Design  of  Flapping-Wing  Micro  Air  Vehicles,” 
Aeronautical  Journal,  Vol.  109,  No.  1098,  Aug.  2005,  pp.  385-393. 

[7]  Chung,  S.  J.,  Stoner,  J.  R.,  and  Dorothy,  M.,  “Neurobiologically 
Inspired  Control  of  Engineered  Flapping  Flight,”  AIAA  Paper  2009- 
1929,  April  2009. 

[8]  Sun,  M.,  and  Xiong,  Y.,  “Dynamic  Flight  Stability  of  a  Hovering 
Bumblebee,”  Journal  of  Experimental  Biology,  Vol.  208,  No.  3,  2005, 
pp.  447—459. 

doi:  10. 1242/jeb. 01407 

[9]  Ellington,  C.  R.  “The  Novel  Aerodynamics  of  Insect  Flight: 
Applications  to  Micro  Air  Vehicles,”  Journal  of  Experimental  Biology, 
Vol.  202,  No.  23,  1999,  pp.  3439-3448. 

[10]  Taylor,  G.  K.,  Bomphrey,  R.  J.,  and 't  Hoen,  J.,  “Insect  Flight  Dynamics 
and  Control,”  AIAA  Paper  2006-32,  Jan.  2006. 

[11]  Died,  J.  M.,  and  Garcia,  E.,  “Stability  in  Omithopter  Longitudinal 
Flight  Dynamics,”  Journal  of  Guidance,  Control,  and  Dynamics, 
Vol.  31,  No.  4,  2008,  pp.  1157-1163. 

doi:  10.25 14/1 .3356 1 

[12]  Doman,  D.  B.,  Oppenheimer,  M.  W.,  and  Bolender,  M.  A.,  "Altitude 
Control  of  a  Single  Degree  of  Freedom  Flapping  Wing  Micro  Air 
Vehicle,”  AIAA  Paper  2009-6159,  Aug.  2009. 

[13]  Wood,  R.  J.,  “The  First  Takeoff  of  a  Biologically  Inspired  At-Scale 
Robotic  Insect,”  IEEE  Transactions  on  Robotics  and  Automation, 
Vol.  24,  No.  2,  2008,  pp.  341-347. 

[14]  Ellington,  C.  P,  “The  Aerodynamics  of  Hovering  Insect  Flight  3: 
Kinematics,”  Philosophical  Transactions  of  the  Royal  Society  of 
London,  Series  B:  Biological  Sciences,  Vol.  305,  No.  1122,  1984, 
pp.  41-78. 

doi:  10. 109 8/rstb.  1984.0051 

[15]  Wood,  R.  J.,  Steltz,  E.,  and  Fearing,  R.  S.,  “Optimal  Energy  Density 
Piezoelectric  Bending  Actuators,”  Sensors  and  Actuators,  A:  Physical, 
Vol.  1 19,  No.  2,  2005,  pp.  476-488. 

doi:  10. 1016/j.sna.2004. 10.024 

[16]  Willmott,  A.  P.,  and  Ellington,  C.  P.,  "The  Mechanisms  of  Flight  in  the 
Hawkmoth  Manduca  Sexta,  1 :  Kinematics  of  Hovering  and  Forward 
Flight,”  Journal  of  Experimental  Biology,  Vol.  200,  No.  21.  1997, 
pp.  2705-2722. 

[17]  Ennos,  A.  R.,  “The  Kinematics  and  Aerodynamics  of  the  Free  Flight  of 
Some  Diptera,”  Journal  of  Experimental  Biology,  Vol.  142, No.  1, 1989, 
pp.  49-85. 

[18]  Doman,  D.  B.,  Oppenheimer,  M.  W.,  and  Sigthorsson.  D.  O., 
“Dynamics  and  Control  of  a  Minimally  Actuated  Biomimetic 
Vehicle,  Part  1:  Aerodynamic  Model,”  AIAA  Paper  2009-6160, 
Aug.  2009. 

[19]  Oppenheimer,  M.  W.,  Doman,  D.  B.,  and  Sigthorsson,  D.  O., 
“Dynamics  and  Control  of  a  Minimally  Actuated  Biomimetic  Vehicle, 
Part  2:  Control,”  AIAA  Paper  2009-6161.  Aug.  2009. 

[20]  Sane,  S.  P.,  and  Dickenson,  M.  H.,  “The  Control  of  Flight  Force 
by  a  Flapping  Wing:  Lift  and  Drag  Force  Production,”  Journal  of 
Experimental  Biology,  Vol.  204,  No.  15,  2001,  pp.  2607-2626. 

[21]  Gradshteyn,  I.  S.,  and  Ryzhik,  I.  M.,  Table  of  Integrals,  Series 
and  Products,  6th  ed..  Academic  Press,  New  York,  2000,  pp.  414- 
415. 

[22]  Doman,  D.,  and  Ngo,  A.,  “Dynamic  Inversion-Based  Adaptive/ 
Reconfigurable  Control  of  the  X-33  on  Ascent,”  Journal  of  Guidance, 
Control,  and  Dynamics,  Vol.  25,  No.  2,  2002,  pp.  275-284. 
doi:10.2514/2.4879 

[23]  Etkin,  B.,  Dynamics  of  Atmospheric  Flight,  Wiley,  New  York,  1972, 
p.  117. 

[24]  Schierman,  J.,  Ward,  D.,  Hull,  J.,  Gandhi,  N.,  Oppenheimer,  M.,  and 
Doman,  D.,  “Integrated  Adaptive  Guidance  and  Control  for  Re-Entry 
Vehicles  with  Flight  Test  Results,”  Journal  of  Guidance,  Control,  and 
Dynamics,  Vol.  27,  No.  6,  2004,  pp.  975-988. 

doi:10.2514/l. 10344 


Approved  for  public  release;  distribution  unlimited. 


16 


